The following work is preliminary and subject to ongoing investigation.
These analyses do not necessarily reflect the point of view of ICCAT or the FAO and in no way anticipate ICCAT or FAO future policy.
Develop a flexible and powerful indicator of stock status for tunas, billfish and sharks that can be tested theoretically and empirically.
Title: ‘Simulation Testing Ecosystem Indicators: Support to ICCAT’s Ecosystem Approach to Fisheries Management (EAFM)’
| Term | Dec 2023 - June 2026 |
| Funding body | International Commission for the Conservation of Atlantic Tunas (ICCAT) |
| Funding stream | FAO Global Environmental Facility (GEF) Common Oceans ABNJ Tuna Project Fund |
| Call No. | # 11903 / 2023 |
| Project Partners | Blue Matter Science Ltd. |
| Blue Matter Team | Tom Carruthers, Adrian Hordyk, Quang Huynh |
| ICCAT Collaborators | Nathan Taylor |
To meet the requirements of the precautionary approach and the ecosystem approach to fisheries management (EBFM), indicators of stock status are needed for secondary species, defined here as those lack sufficient data or capacity to conduct routine stock assessments (see Carruthers et al. 2024 for a detailed problem statement). These indicators must be theoretically sound and be validated empirically.
The EcoTest project aims to use simulation modelling to identify data and algorithms that can inform stock status of secondary species and then validate these indicators empirically in cases where there have been defensible stock assessments.
A well-documented, defensible and transparent framework is needed to support tactical decision making to move beyond the single species assessment paradigm and make progress towards the essential goals of EBFM.
Previous work synthesized the dynamics of six stock assessments and consolidated those into a multi-species, multi-fleet operating model (Huynh et al. 2022). Those operating models were used to simulate various future scenarios for fishing and stock dynamics to develop Indicator 1: an indicator system specific to the longline case study that could predict stock status for the 6 species without undergoing a full stock assessment (Carruthers et al 2024).
The key advance of this earlier work was to use deep learning to train neural networks on the posterior-predicted data of a large number of simulations (a training dataset). The ETest package is built around an evoluation of that approach, using simulated posterior predicted data from a much more general set of simulations that include population and fishing dynamics for a wide range of pelagic highly migratory species that are caught in Atlantic fisheries.
These indicators focus on the prediction of stock status relative to a productive level expressed as spawning stock biomass relative to spawning stock biomass corresponding to MSY levels. This is a key reference level for evaluating stock status for the majority of stock assessed by ICCAT. Expert working groups also confirmed that there is interest in developing indicators that can accommodate various numbers of fleets and stocks, and that these should also apply to varying data availability scenarios with respect to length, age and relative abundance data.
EcoTest can accept a range of data inputs (features) that may be available for secondary species derived from data streams such as total annual catches, CPUE indices, mean length in catches and mean age in catches. These can be provide along with life-history (e.g. somatic growth, natural survival) and fleet-specific (e.g. selectivity) characteristics to potentially improve reliability in interpreting these data inputs.
Where multiple fleets can be identified, data inputs can be provided for up to 3 individual fleets.
Where multiple stocks can be assumed to have been exploited by a common fishery with comparable exploitation history, multi-stock data inputs can be used including catch ratios of one stock to another.
A training dataset was generated from a large number of multi-stock, multi-fleet simulations spanning a range of life-histories and fishery dynamics typical for tuna, billfish and sharks. In each simulation more than 400 data inputs were calculated along with the known, true simulated stock status (spawning biomass relative to spawning biomass at MSY levels). The data inputs relate to trajectories, levels and variability of data such as catches, CPUE indices and mean length in catches and are processed and formatted in such a way that they can be calculated for any fishery (e.g. mean lengthin catches is phrased as a percentage of asumptotic length, slopes in data are slopes in log data etc).
Some data inputs relate to data from fishing fleets, some relate to fishery characteristics, and others relate to the life hsitory of the stock(s) (see Table 1).
The large training dataset is subsetted to include only those data inputs provided by the user. An artificial neural network is then trained on the substted training dataset.
During training a validation dataset is used to check for model overparameterization and confirm predictive ability for independent data not used by the training algorithm.
After training, the expected indicator performance is evaluated according to a completely independent testing dataset.
This step can be seen as (a) a theoretical check that there is information in the features (input data) provided to estimate stock status and (b) a quantitative theoretical evaluation of the expected performance of the indicator.
A result of this step is a fitted model and also a standardized version of the user’s dataset that has the same standardization process as the training dataset.
The product of step 3 (fitted model and standardized version of the dataset) can be used to estimate stock status.
If a historical trend in empirically estimated stock status is provided, the estimates of the indicator can be compared with other empirical estimates such as those arising from stock assessments.
The Keras, Tensorflow and Miniconda libraries on which EcoTest methods are based are highly sensitive to the version of R and RStudio.
It is therefore highly recommended that users install the very latest version of R and RStudio before continuing (you will almost certainly have difficulties otherwise).
R can be downloaded from the R Project for Statistical Computing webpage.
Currently the ETest package is compatible with R version R version 4.5.1 (2025-06-13 ucrt).
RStudio can be downloaded from the Posit webpage.
Currently the ETest package is working in RStudio version 2025.05.1 Build 513.
The ETest package has only been tested on a Windows x64 system.
Open RStudio and run the following lines of code at the R command prompt:
# Keras, Tensorflow, Miniconda
install.packages("keras3")
library(keras3)
install_keras()
# EcoTest package ETest
install.packages("remotes")
remotes::install_github('blue-matter/ETest')
To test your installation run the following code
library(ETest)
test = train_ind(Shark_1_data)
If your packages are working you should see something like this:
Lets just assume we have some data on catches, recent CPUE and a bit of life history info. Remember, these names have to match the EcoTest codes for indicator features (see Table 1 below).
library(ETest)
dat = data.frame(
K_s1 = 0.25, # von B. somatic growth (stock 1)
M_K_s1 = 1.05, # Instantaneous natural Mort. / von B. K (stock 1)
L50_Linf_s1 = 0.74, # Length at 50% maturity (spawning fraction) / asymptotic length (L-infinity) (stock 1)
maxa_s1 = 11, # Age at 5% survival (stock 1)
L5_L50_s1_f1 = 0.2, # Shortest length at 5% selectivity relative to length at 50% maturity (stock 1)
C_rel_s1_f1 = 0.67, # Current catches relative to historical average (stock 1, fleet 1)
C_g10_s1_f1 = 0.01, # Slope in log catches over the last 10 years (stock 1, fleet 1)
I_rel_s1_f1 = 0.43, # Current CPUE relative to historical average (stock 1, fleet 1)
I_g20_s1_f1 = -0.09, # Slope in log CPUE over the last 20 years (stock 1, fleet 1)
ML_rel_s1_f1 = 0.9 # Current mean length in catches relative to historical average (stock 1, fleet 1)
)
input = list(dat=dat)
my_indicator = train_ind(input)
You can see that the performance of this indicator isn’t stellar however it correctly classifies stocks below half SSBMSY (an ICCAT limit reference point) and above SSBMSY (an ICCAT target reference point) in about 60% of simulations.
For now we are going to assume we find this acceptable (you might well not!).
We can now obtain a point estimate of stock status using our indicator and the dataset provided:
my_pred = pred_ind(my_indicator)
my_pred
Since only one set of data were provided there is only a point estimate of stock status. However most datasets can be sampled to provide stochastic inputs (bootstrapping etc).
For the purposes of demonstration lets just create some samples of data so you can plot indicator estimates that include uncertainty:
nsamp = 50 # 50 stochastic draws of the dataset
nobs = ncol(dat) # data set has nobs features (data types)
log_norm_err = rlnorm(nobs*nsamp,0,0.2) # demo error
data_repped = rep(as.numeric(dat),each=nsamp) # duplicate data for nsamp rows
dat_s = array(data_repped * log_norm_err,c(nsamp,nobs)) # combine data with errors
dat_s = as.data.frame(dat_s) # make into data frame
names(dat_s) = names(dat) # make sure to copy over labels
input_s = list(dat=dat_s) # make an input object where position one in the list is the data
my_indicator_s = train_ind(input_s) # train the indicator
my_pred_s = pred_ind(my_indicator_s) # make a status prediction for each of the 50 samples of data
hist(my_pred_s, xlab="Estimated SSB/SSBMSY", main="Demonstration Stochastic Indicator", col="cornflowerblue", border="white")
In this example we will process some example time series data, train an indicator and estimate stock status using it.
The ETest package comes with an example dataset ‘Some_data’ for these worked examples. The dataset has catches, CPUE indices, mean length, life history and selectivity information for three fleets simultaneously fishing three fleets. We are not necessarily going to assume we have all these data, but they are there to experiment with.
To make creating time series data inputs easier, ETest has a function ts_features() which extracts the levels and slopes processed in the same way as the training dataset.
stock = 1 # Lets just try this for one stock and ...
fleet = 1 # ... one fleet.
Catch = Some_data$Catch[stock,fleet,] # Catch data is a 3D array [stock, fleet, year]
plot(as.numeric(names(Catch)),Catch) # what do these catch data look like?
Catch_f = ts_features(Catch, "C") # lets capture the rel, g5, g10, g20 and g40 features
Catch_f
Index = Some_data$Index[stock,fleet,] # CPUE index data is a 3D array [stock, fleet, year]
Index_f = ts_features(Index, "I") # captures the rel, g5, g10, g20 and g40 features
ML = Some_data$Mean_length[stock,fleet,] # Mean length data
ML_f = ts_features(ML,"ML") # captures the rel, g5, g10, g20 and g40 features
fleet_data = as.data.frame(c(Catch_f, Index_f, ML_f)) # combine fleet data
names(fleet_data) = paste0(names(fleet_data),"_s",stock,"_f",fleet) # label correctly
LH_sel_data = data.frame(K_s1 = Some_data$K[stock], # combine life history and selectivity info
M_K_s1 = Some_data$M[stock]/Some_data$K[stock], # ratio of M to K
maxa_s1 = Some_data$maxa[stock], # age at 5% natural survival
L5_L50_s1_f1 = Some_data$L5[stock,fleet]/Some_data$L50[stock], # shortest length at 5% selectivity relative to L50
LFS_L50_s1_f1 = Some_data$LFS[stock,fleet]/Some_data$L50[stock]) # length at 95% selectivity relative to L50
input2 = list(data=cbind(LH_sel_data, fleet_data)) # Combine inputs into a single input object
Now that we have our dataset, extracted from three time series (Catch, a CPUE index and Mean Length) in addition to selectivity and life history parameters, we can train and use the indicator:
my_indicator2 = train_ind(input2) # train indicator, standardize input data
And then get an estimate of stock status:
pred_ind(my_indicator2) # estimate stock status from trained indicator and submitted dataset
The ETest package comes with a number of real fishery datasets loaded. You can list these using the R function data():
data(package="ETest")
There is an object TD which is the full training dataset. You can always look at that if you want to check the magnitude and range of data inputs used for training the neural networks.
The other objects are labelled _io, _data and _retro.
_io objects such as Shark_1_io are list objects that contain all of the inputs and outputs to a Stock Synthesis 3 assessment.
_data objects are list objects that contain ETest formatted features obtained from those _io stock synthesis runs (position 1, ‘inputs’) and also the stock synthesis estimated stock status (position 2, ‘Brel’).
_retro objects are a list of _data objects npeels long which strip back the data from the stock synthesis runs to allow retrospective behavior of indicators to be evaluated.
Using these objects it is simple to conduct indicator training and prediction, including retrospective analysis:
# Single stock status estimate using data to most recent year:
Shark_1_ind = train_ind(Shark_1_data)
Shark_1_pred = pred_ind(Shark_1_ind)
hist(Shark_1_pred)
# Retropsective stock status estimate using npeels of the SS3 data:
Shark_1_ind_retro = train_ind(Shark_1_retro)
Shark_1_pred_retro = pred_ind(Shark_1_ind_retro)
To obtain all relevant outputs for a current Stock Synthesis model you will a directory where an assessment has been run and you will need to install the latest version of the r4ss package.
wm_io = all_ss3("C:/myassessments/ATL_white_marlin_2025")
Once you have the assessment ripped into R, you need to pick what fleets and indices you wish to extract. You can use the ss_names() function to list the fleets / surveys. You then have to select the relevant fleets and corresponding cpue indices and then make your formatted EcoTest indicator features using SS_2_ET():
ss_names(wm_io)
wm_data = SS_2_ET(wm_io, Fnam = c("JPN_LL_N","US_PL_N"), Inam = c("Surv_JPN_LL_N","Surv_US_PL_N"))
The wm_data object now has features for three fleet data sets: the Japanese longline north the US pole and line north and all other fleets combined. Selectivities for all other fleets combined are weighted by historical average catches.
The function SS_2_ET() samples from time series and parameters to provide a stochastic feature dataset in order to obtain variance in stock status estimates.
You can now run the indicator for that dataset:
wm_ind = train_ind(wm_data)
wm_pred = pred_ind(wm_ind)
or alternatively you can do the same steps but conduct a retrospective analysis using the function SS_2_ET_Retro():
wm_retro = SS_2_ET_Retro(wm_io, Fnam = c("JPN_LL_N","US_PL_N"), Inam = c("Surv_JPN_LL_N","Surv_US_PL_N"),npeels=8)
wm_ind_retro = train_ind(wm_retro)
wm_pred_retro = pred_ind(wm_ind_retro)
In general, fisheries stock assessment has failed to incorporate information across multiple species in the estimation of stock status. This is a missed opportunity.
For example if the denominator of CPUE indcies of abundance (effort) is comparable in trend among two species, their catch ratio informs their relative depletion (e.g. if the catch ratio of stock 1 : stock 2 has declined by 50%, then this infers that stock 1 is at half the vulnerable stock size as stock 2 (assuming you standardized your data correctly).
This means that depletion on one stock may be informed by data informing depletion on the other stock. As the number of stocks increases (that adhere to this assumption of approximately comparable effort trends) an increasibly narrow set of depletion conditions may prevail. To demonstrate this lets see how well only catch data can inform stock status across three species:
fleet = 1 # ... one fleet.
Catch_f1 = ts_features(Some_data$Catch[1,fleet,],"C")
Catch_f2 = ts_features(Some_data$Catch[2,fleet,],"C")
Catch_f3 = ts_features(Some_data$Catch[3,fleet,],"C")
cvec = c(unlist(Catch_f1),unlist(Catch_f2),unlist(Catch_f3))
msdat = data.frame(matrix(cvec,1))
names(msdat) = paste0(names(cvec),"_s",rep(1:3,each=5),"_f1")
input = list(dat=msdat)
catch_ind = train_ind(input)
Since there is nothing providing an indication of absolute depletion the indicator based on catches alone performs very poorly.
fleet = 1 # ... one fleet.
Index_s1 = ts_features(Some_data$Index[1,fleet,],"I")
Index_s2 = ts_features(Some_data$Index[2,fleet,],"I")
Index_s3 = ts_features(Some_data$Index[3,fleet,],"I")
Index_rels = data.frame(I_rel_s1_f1 = Index_s1$I_rel,
I_rel_s2_f1 = Index_s2$I_rel,
I_rel_s3_f1 = Index_s3$I_rel)
wIdat = list(dat = cbind(msdat, Index_rels))
catch_wI_ind = train_ind(wIdat)
Adding just the current index levels relative to historical means, tidied the indicator up somewhat.
Now lets find out how much of that was just the index levels:
catch_onlyI_ind = train_ind(list(dat=Index_rels))
As expected, the relative index levels were doing the majority (but not all) of the lifting.
< coming soon!>
< coming soon!>
< coming soon!>
Table 1. Indicator Features (data types)
| Feature | Type | Description |
|---|---|---|
| K | Stock | von Bertalanffy growth parameter (annual) |
| M_K | Stock | The ratio of natural mortality rate to von B. K |
| maxa | Stock | Maximum age. The age corresponding to 5% natural survival (the age before which 95% of individuals have died from natural mortality) |
| L50_Linf | Stock | The ratio of length at 50% maturity (spawning fraction) to asymptotic length (von B. L-infinity) |
| L5 | Fleet | The shortest length at 5% selectivity phrased as a % of asymptotic length (L-infinity) |
| LFS | Fleet | The shortest length at 95% selectivity phrased as a % of asymptotic length (L-infinity) |
| VML | Fleet | The vulnerability at maximum length (‘dome-shapedness’) |
| I_rel | Data | The ratio of current index level to mean historical index level |
| I_g5 | Data | The slope in log index over the most recent 5 years |
| I_g10 | Data | The slope in log index over the most recent 10 years |
| I_g20 | Data | The slope in log index over the most recent 20 years |
| Isd1 | Data | Mean absolute residual error in a log CPUE index from a loess smoother with effective number of parameters of 40% index length |
| Isd2 | Data | As Isd1 but ENP of 20% |
| Isd3 | Data | As Isd1 but ENP of 10% |
| C_rel | Data | The ratio of current catches relative to mean historical levels |
| C_g5 | Data | The slope in log catches over the most recent 5 years |
| C_g10 | Data | The slope in log catches over the most recent 10 years |
| C_g20 | Data | The slope in log catches over the most recent 20 years |
| Csd1 | Data | Standard deviation in residual error around log catches from a loess smoother with ENP 20% catch series length |
| CF | Data | The fraction of total historical catches (taken by this fleet). |
| ML_rel | Data | The ratio of current mean length in catches relative to historical average mean length in catches. Mean length is phrased as a % of asymptotic length (L-infinity) |
| ML_g5 | Data | The slope in log of mean length of catches over the most recent 5 years |
| ML_g10 | Data | The slope in log of mean length of catches over the most recent 10 years |
| ML_g20 | Data | The slope in log mean length of catches over the most recent 20 years |
| ML_Linf | Data | The ratio of current mean length in catches relative to asymptotic length (von B. L-infinity) |
| ML_L50 | Data | The ratio of current mean length in cathes to length at 50% maturity (spawning fraction) |
| MA_rel | Data | The ratio of current mean age in catch relative to mean historical levels |
| MA_g5 | Data | The slope in log mean age in catch over the most recent 5 years |
| MA_g10 | Data | The slope in log mean age in catch over the most recent 10 years |
| MA_g20 | Data | The slope in log mean age in catch over the most recent 20 years |
| MV_rel | Data | The standard deviation in observed caught lengths in the current year relative to the historical average. |
| MV_g5 | Data | The slope in log standard deviation in catches over the most recent 5 years |
| MV_g10 | Data | The slope in log standard deviation in catches over the most recent 10 years |
| MV_g20 | Data | The slope in log standard deviation in catches over the most recent 20 years |
| FM_rel | Data | The ratio of current fraction mature fish in catches relative to historical mean fraction of mature fish in catches |
| MV_g5 | Data | The slope in log of fraction mature fish in catches over the most recent 5 years |
| MV_g10 | Data | The slope in log of fraction mature fish in catches over the most recent 10 years |
| MV_g20 | Data | The slope in log fraction mature fish in catches over the most recent 20 years |
| CR | Data | The current catch ratio (first stock divided by second stock). R(y) = C(1,y) / C(2,y). CR is R from the most recent year, R(ny) |
| CR_mu | Data | The catch ratio (first stock divided by second stock) over the first 20 years R[1:20] |
| CR_rel | Data | CR / CR_mu |
| CR_s5 | Data | The slope in log CR over the most recent 5 years |
| CR_s10 | Data | The slope in log CR over the most recent 10 years |
| CC20 | Data | Correlation of the in residual fit of a loess smoother to log catches with ENP = 20% among stocks over the most recent 20 years. |
| CC40 | Data | As CC20 but over historical years between 21 and 40 years ago. |
| CC60 | Data | As CC20 but over historical years between 41 and 60 years ago. |
The openMSE website contains extensive documentation about openMSE and approaches such as RCM.
Zhou’s intro to ANN and machine learning
Tensorflow: Platform for machine learning
Nvidia CUDA Toolkit: GPU accellerated applications
Carruthers, T.R. 2018. A multispecies catch-ratio estimator of relative stock depletion. Fisheries Research. 197: 25-33, https://doi.org/10.1016/j.fishres.2017.09.017.
Huynh, Q., 2025. Rapid Conditioning model. Available from https://openmse.com/tutorial-rcm/
Juan-Jordá, M.J., Murua, H., Arrizabalaga, H., Dulvy, N.K., and Restrepo, V. 2018. Report card on ecosystem-based fisheries management in tuna regional fisheries management organizations. Fish Fish. 19: 321-339.
EcoTest is funded by ICCAT as part of a commitment to the Global Environmental Facility (GEF) Common Oceans ABNJ Tuna Project Fund.
Many thanks to The Ocean Foundation for funding Phases I and II of the EcoTest framework.
Many thanks to the ICCAT Ecosystem working group for comments on earlier versions of EcoTest: Nathan Taylor, Alex Hanke, Sachiko Tsuji, Guillermo Diaz, Diego Alvarez, Laurie Kell, Maria Juan Jorda, Andres Domingo.
Blue Matter: Tom Carruthers, Adrian Hordyk, Quang Huynh
OpenMSE was developed with support from the Natural Resources Defense Council (NRDC), the Gordon and Betty Moore Foundation, the Packard Foundation, the Marine Stewardship Council, Fisheries and Oceans Canada (DFO), the U.S. National Oceanic and Atmospheric Administration, the International Commission for the Conservation of Atlantic Tunas (ICCAT) and The Ocean Foundation.
OpenMSE continues to be developed with the support of the The Ocean Foundation.